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Within the established theoretical framework of quantum mechanics, interference always occurs 
between pairs of trajectories. Higher order interferences with multiple constituents are, however, 
excluded by Born’s rule and can only exist in generalized probabilistic theories. Thus, high-precision 
experiments searching for such higher order interferences are a powerful method to distinguish 
between quantum mechanics and more general theories. Here, we perform such a test in optical 
multi-path interferometers. Our results rule out the existence of higher order interference terms 
to an extent which is more than four orders of magnitude smaller than the expected pairwise 
interference, refining previous bounds by two orders of magnitude. This establishes the hitherto 
tightest constraints on generalized interference theories. 


Since arising almost a century ago, quantum mechan¬ 
ics has long become an established paradigm for the de¬ 
scription of nature on a submicroscopic scale. It is at the 
basis of an enormous variety of present and potential fu¬ 
ture applications, such as quantum communication [1, 2], 
quantum computation [3-5] and protocols like entangle¬ 
ment swapping [6] or teleportation [7]. However, all these 
applications rely ultimately on interference and entangle¬ 
ment, which can be alternatively explained by theories 
sharing only some fundamental features with quantum 
mechanics, such as the superposition principle or proba¬ 
bilistic predictions for outcomes, and yet differing from it 
in other aspects. In order to distinguish between quan¬ 
tum theory and such alternatives one needs to design ded¬ 
icated experiments. The situation may be compared with 
the time before the first Bell test experiments had been 
performed. Until then one could explain all quantum me¬ 
chanical phenomena with a local hidden variable theory. 
It was required to first state Bell’s theorem [8] and then 
to perform dedicated experiments with space-like sepa¬ 
rated laboratories to exclude the alternative. Only last 
year all experimental loopholes were finally closed [9-11]. 
Another example is the experiment that distinguishes be¬ 
tween quaternion and complex (standard) quantum the¬ 
ory [12, 13]. Here we focus on an experimental test ca¬ 
pable of discerning between quantum mechanics and its 
generalizations exhibiting higher order interference. 

The probabilistic nature of quantum theory is stated 
by Born’s rule [14], i.e. that the probability density 
P(r, t) for an observation of a quantum object at a certain 
time t and a certain position r is given by the absolute 
square of its wavefunction 'I'(r, f): 

P(r,<) = T*(r,f)4'(r,t) = |«'(r,f)l^ . (1) 

As a consequence of Born’s rule and quantum superposi¬ 
tion, interference can take place even for single particles 
[15]. For concreteness, consider an interferometer with 
multiple non-overlapping paths k = A, B,C,... which 
superpose in some output port to the hnal wavefunction 


4' = Y.k'^k- Eq. (1) implies: 

P<A,t) = X + (2) 

k k<.l 

with pairwise (first-order) interference terms lu = 
'kfeT'* -I- C.C., depending on the relative phase between 
the two paths k and 1. Thus, one obtains interference 
terms that always originate from pairings of paths, but 
no higher order interferences involving more than two 
paths at once. 

In this vein, one can use the presence or absence of 
higher-order interferences as an experimental probe of 
the current framework of quantum mechanics. First de¬ 
veloped by Sorkin in the context of a measure theory 
on spacetime [16], one can define a hierarchy of interfer¬ 
ence terms. In a 3-path interferometer with individually 
blockable paths A, B, C, where Pabc is the probability 
to Hnd a particle in the output port of the interferometer 
if all paths are open, Pab for only paths A and B being 
open, etc. The so-called second-order interference term 

I ABC = Pabc — Pab — Pac — Pbc + Pa + Pb + Pc (3) 

should be zero, independent of the individual phases and 
powers in each interferometer arm, due to Eq. 2. Con¬ 
versely, a signihcant deviation from Iabc = 0 would 
indicate the existence of higher-order interferences and 
contradict conventional quantum theory. Note that the 
definition (3) accounts for deviations from the standard 
theory in a model-independent way. 

In any experiment with discrete particles, the proba¬ 
bility P will be proportional to the detected particle flux 
p. Therefore a directly measurable quantity 

€3 = PABC - PAB - PAC “ PBC +PA+PB+PC -PO (4) 

can be defined [17]. In this expression, for example is 
the detected particle flux at the output when only paths 
A and B are open. The background term po gives the 
measured signal when all paths are blocked, account¬ 
ing for detector dark current/dark counts. For better 
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comparison of the results with the expected behavior, 
one can introduce the normalized quantity K 3 = 
measuring the ratio of hypothetical second-order interfer¬ 
ence to the sum of the expected first-order interference, 
1^3 = I Iab I + 1 1 AC I + 1 I ■ For interferometers with more 

than three paths, the higher (third, fourth,...)-order in¬ 
terference terms (k4, K5,...), which of course are also zero 
in standard quantum theory, can be defined accordingly. 

Different experiments have been realized previously to 
obtain an upper bound on the modulus of the second- 
order interference term. These experiments were imple¬ 
mented in optics [1, 18, 19] as well as via nuclear mag¬ 
netic resonance (NMR) in molecules [21], delivering re¬ 
sults which were all in accordance with the expectation 
K 3 = 0. The NMR experiment provided the hitherto 
tightest constraint with K 3 = 0.001 ± 0.003. 

As is the case for any such null-test experiment, the 
tightness of the bound and, thereby, the strength of any 
conclusions to be drawn about the foundations of the the¬ 
ory depend on the measurement uncertainties. In pre¬ 
vious optical 3-path interferometers, the precision was 
mostly limited by the phase stability of the interferom¬ 
eter, while the accuracy suffered from detector nonlin¬ 
earities [1]. In this work, we present a greatly improved 
multi-path experiment, namely a stabilized 5-path inter¬ 
ferometer with single photons, with which we are not only 
able to tighten the bound on second-order interference 
by two orders of magnitude, but also measure third and 
fourth-order interference terms. The 5-path interferome¬ 
ter has the additional advantage of permitting the acqui¬ 
sition of more statistics for the second- and third-order 
interference term since it consists of ten 3-path interfer¬ 
ometers and five 4-path interferometers. The systematic 
error of detector nonlinearities is taken into account by 
separate detector calibration and full quantum state to¬ 
mography of the produced 5-dimensional qudit state. 

A schematic drawing of our setup can be seen in Fig. 1. 
We employed light sources in three different regimes: 
classical (C), semi-classical (SC) and quantum (Q). In 
the classical regime we used a continuous-wave single fre¬ 
quency laser at 808 nm power-stabilized to 1 mW with 
relative fluctuations smaller than 0 . 1 % over the com¬ 
plete measurement time of several days by using a liquid 
crystal noise eater (Thorlabs LCC3112). We used a sin¬ 
gle photon source to perform measurements in the other 
regimes. Photon pairs at 808nm are produced via type- 
II spontaneous parametric down conversion (SPDC) in 
a 10 mm long periodically poled potassium titanyl phos¬ 
phate (ppKTP) crystal which is pumped by a blue laser 
(404nm). The orthogonally polarized photons are sep¬ 
arated on a polarizing beam splitter (PBS). We collect 
6 X 10® single photons per second in each of the outputs 
in single mode fibers and we get 10 ® pairs per second 
at 4mW pump power. One of the photons serves as 
a heralding photon, whereas the other is sent through 
our multi-path interferometer. Therefore we have two 



FIG. 1. Experimental setup. The light source is either given 
by a power stabilized laser or single photons (both at 808 nm) 
produced via SPDC in a ppKTP crystal pumped by a blue 
laser. The interferometer consists of two diffracitve beam¬ 
splitters and two lenses. Shutters and phase plates in each of 
the paths allow independent manipulation. The inset shows 
the dimensions and separations of the shutters and the 1/e^ 
intensity diameters of the beams. 


possibilities to conduct the measurement with single 
photons: either free running, where all photons trans¬ 
mitted through the interferometer are counted (yield¬ 
ing a thermal photon number distribution) in the semi- 
classical regime or conditioned, where only photons are 
counted if there is a heralding photon (producing a sub- 
Poissonian distribution) [22] in the quantum regime. All 
light sources were linearly polarized. 

The interferometer is a Mach-Zehnder 5-path interfer¬ 
ometer consisting of a diffractive beam splitter (a diffrac¬ 
tive optical element - Holoeye DE 263 - modulating 
the incident light via a micro-relief surface) which cre¬ 
ates five almost equally powerful beams, collimated by 
a lens (/ = 150 mm). A shutter assembly serves to 
block or unblock each of the five beams individually, 
phase plates (glass plates with a thickness of 0.15 mm 
and anti-reflection coated for 808 nm) mounted on mo¬ 
torized rotation stages in all of the five beams allow us 
to set the phase of each path independently. The abso¬ 
lute angular repeatability is 0.005°, which corresponds 
to 7 r /1000 in phase. An identical pair of lens and diffrac¬ 
tive element sends the resulting beam onto a detector. 
The interferometer is designed in 4/-configuration and 
the individual beams are separated by 5 mm, therefore 
the overall dimensions are (60 x 2) cm^. For detect¬ 
ing single photons we used SPCM-AQRH-12-FC single 
photon counting modules from Perkin Elmer followed by 
a quTAU time-to-digital converter from qutools GmbH. 
This system has a deadtime of (33.85 ± 0.31) ns and 
(150 ± 18) dark counts per second. The laser radiation 
is detected by a Physimetron A139-001 photoreceiver 
based on a Si-photodiode (Hamamatsu S2386-18K) and 
a 1 MV /A transimpedance ampliher, read out by an Ag- 
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FIG. 2. Mean values of the measured powers for the different 
path-combinations with classical light. Extreme outliers have 
already been removed. 


ilent 34410A multimeter. This detection system has a 
low maximum nonlinearity of less than 35 ppm [3]. One 
measurement set consists of the 2® = 32 different pos¬ 
sible open/close combinations of the five paths. These 
combinations were measured in random order to reduce 
the influence of any memory effects of the detector and 
of drifts of the source. To obtain data with comprehen¬ 
sive statistics we recorded several thousand measurement 
sets within a total measurement time of several days. 
The whole interferometer is shielded against air motion 
and stray light as well as passively and actively tempera¬ 
ture stabilized with a PI controller (Wavelength Electron¬ 
ics HTC1500) and heating mats to a root-mean-square 
fluctuation < 0.02 K /24 h (The temperature was mon¬ 
itored with a PTIOOO resistance thermometer). Addi¬ 
tionally, the phases are actively stabilized by optimizing 
the phase-plate position after 100 measurement cycles 
towards maximally constructive interference of all two- 
path combinations. This point in phase space were cho¬ 
sen for convenience of alignment and because small phase 
changes lead only in second order to deviations in output 
power. This results in good phase stability over the whole 
measurement time. By comparing the fluctuations and 
drifts of the single-path power with the multi-path pow¬ 
ers (see Fig. 4 and Fig. 5 of the Supplemental Material 
[24]), which have the same order of magnitude, we found 
that phase uncertainty plays a minor role compared to 
power noise from the light source. 

The resulting average powers of the different path com¬ 
binations can be seen in Fig. 2 for the measurement with 
the power stabilized laser [24]. We Altered the data for 
extreme outliers (resulting from shutter failure) accord¬ 
ing to Grubbs’ test for outliers (with a significance level 
of 99%) [25, 26]. After that the largest relative stan¬ 
dard deviation of the various classical signals is 0.3% 
for 5618 measurement sets recorded within 68 h. For 
the semi-classical (quantum) single photon measurement 
with 1912 measurement sets the largest standard devia¬ 
tion was measured to be 3.6% (15.5%) over a measure¬ 


FIG. 3. Histogram of £ 3 /( 53 ) values in the classical regime for 
the 3-path subset {A,B,C}. The blue line is a Gaussian fit to 
the distribution of the data and the mean value is indicated 
by the vertical black line, whose width corresponds to the 
standard error of the mean (one standard deviation). 

ment time of 88 h. These higher values result mainly from 
shot noise and from the fact that the power of the blue 
pump laser was not stabilized. 

Due to the anti-correlation between the numerator e 
and the denominator S in the definition of k, a bias to¬ 
wards positive values can arise from random fluctuations 
in the data when calculating k for every shutter cycle 
and averaging over the data sets. However, calculating 
the averages of numerator and denominator in the defi¬ 
nition of K separately, 

J= 3,4,5 (5) 

eliminates their correlations and yields an unbiased es¬ 
timator of the higher-order interference terms. Indeed, 
one can show that error sources, which typically occur 
in interference experiments, such as power fluctuations 
of the photon source, countrate fluctuations of the detec¬ 
tors (Poissonian photon counting uncertainties), detec¬ 
tor/electronic noise, coherent phase fluctuations as well 
as incoherence, have no systematic effect on the measure¬ 
ment outcome [27]. 

For each of the measured 32-tuples we calculate £3,4,5 
and () 3 , 4 , 5 . A histogram plot of the measured ensemble 
of £ 3 /( 53 ) in the classical regime is shown in Fig. 3 [28]. 
After averaging across all possible path combinations one 
obtains the mean values and associated uncertainties pre¬ 
sented in table I. 



(^ 3 ) 

(k4) 

(«5) 

classical (xlO””) 

9.7±0.1 

2.7 ±0.2 

0.3 ±0.3 

semi-classical (xl0“^) 

-9.9 ± 1.8 

-5.1 ±2.1 

-3.8 ±3.9 

qnantum (xl0“^) 

-1.1 ± 1.6 

0.3 ± 1.8 

-2.6 ± 2.9 


TABLE I. Mean values of the measured higher-order interfer¬ 
ences and their standard errors in the classical, semi-classical 
and quantnm regimes. 

It was recently shown that near-field effects in slit- 
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based measurements, where the relevant dimensions are 
just one or two orders of magnitude larger than the wave¬ 
length A, can lead to an apparent higher-order interfer¬ 
ence and, therefore, bias the experiment [29, 30]. How¬ 
ever, in our interferometer these effects are of negligible 
influence, due to the macroscopic separation of the paths 
and beam width, exceeding A by 3 to 4 orders of magni¬ 
tude (see inset in Fig. 1). Instead, the main systematic 
uncertainty in our experimental configuration arises from 
the nonlinearity of the detectors. Real detectors usually 
have a nonlinear response function, which means that the 
recorded value (voltage, photon counts,...) is not linear 
in the incident power or photon flux, but biased differ¬ 
ently for different optical powers. This biases the value 
of K, as we measure light powers varying over more than 
one order of magnitude. Here, the bias arises mainly due 
to nonlinearities in the electronics of our photoreceiver 
and due to deadtime in the single photon detector. To 
take this error into account it is useful to fully character¬ 
ize our 5-path interferometer, which can be described as 
a 5-dimensional qudit state. Therefore, we additionally 
performed complete quantum state tomography [31, 32]. 
The density matrix p was numerically reconstructed from 
single- and two-path measurements with defined phases 
via direct reconstruction. The phases are calibrated via 
scanning the classical two-path laser interference. We 
used the direct reconstruction instead of a maximum like¬ 
lihood estimation to avoid systematic deviations in the 
state reconstruction, which have recently been shown to 
arise due to the constraint of physicality in maximum 
likelihood estimates [33]. The real and imaginary parts 
of the resulting density matrix are shown in Fig. 4. We 
calculated trp^ = 0.74; the deviation from 1 (a pure 
state with no which-path information) can be attributed 
to an imperfect overlap of the five beams at the second 
beamsplitter. While this degree of coherence in the inter¬ 
ferometer must be determined for an accurate prediction 
of the influence of the nonlinearities, its actual value has 
no systematic impact on the Sorkin experiment. The 
effect of the nonlinearity on the reconstruction is neg¬ 
ligible for two reasons: the ratio of the photon fluxes 
for the different measurement settings is much smaller 
than in the measurements contributing to the evaluation 
of K. More importantly, deviations in the density ma¬ 
trix do not produce a systematic effect on the expected 
higher-order interference, as k = 0 holds for all states in 
quantum theory. 

From the density matrix it is possible to calculate the 
expected powers for the different settings of the shut¬ 
ters in the Sorkin experiment. We found good accor¬ 
dance with our measurement data [34], suggesting that 
the tomography produces an accurate description of the 
interferometer. The nonlinearities of both detectors have 
been characterized in separate experiments [3]. Applying 
them to the powers/count rates predicted from the den¬ 
sity matrix yields small corrections of these powers (rel- 


Re (p) Im (p) 



FIG. 4. The density matrix p of the 5-dimensional qudit state 
in onr interferometer, with trp^ = 0.74. 
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FIG. 5. Final result. R gives the difference between the exper¬ 
imentally measnred (k) and the value expected from detector 
nonlinearities Kth for the measurements in the three regimes. 
The order of R increases along the horizontal axis and the 
error bars indicate one standard deviation. 

ative change < 0.03% for the laser powers and < 0.5% 
for the unheralded single photon rates). 

One can then calculate the apparent higher-order in¬ 
terferences Kth, which would be expected in the Sorkin 
measurements, from these corrected data. The differ¬ 
ences between the experimentally measured higher-order 
interferences and the expected values due to the nonlin¬ 
earities K = (k) — Kth give corrected higher-order inter¬ 
ferences as the final results, which can be found in table 
H. 



ka 

k;4 

ks 

classical (xlO”'') 

0.0 ±3.1 

4.3 ±4.4 

4.2 ±5.1 

semi-classical (xl0“^) 

1.3 ± 1.8 

-1.6 ±2.1 

-3.8 ±4.0 

quantum (xl0“^) 

0.0 ± 1.6 

0.6 ± 1.8 

-2.7 ± 2.9 


TABLE II. R = (k) — Kth is the nonlinearity-corrected higher- 
order interference for all measurement regimes. All these val¬ 
ues are within one standard deviation of the expected zero 
value. 

Note that in case of the heralded single photon data, 
we did not calculate an explicit prediction for Kth because 
the nonlinearity model is quite involved in this case. In¬ 
stead we used the model to correct the raw experimental 
data, in order to obtain k [35]. A final summary of all 
the different kj values is presented in Fig. 5. One Ends 
that all these values are within one standard deviation of 
the expected zero value. 

The optical 5-path interferometer presented in this 
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work permitted us to experimentally confine the allowed 
domain of second order interference to an uncertainty of 
3 X 10“® in the classical light regime. This is two orders 
of magnitude tighter than the bounds obtained from the 
most precise experiments in any system to date. The un¬ 
certainties in the semi-classical and quantum regimes of 
2 X 10 ^ and 2 X 10 respectively, are also much lower 
than what has been reported before [1, 17]. This new 
level of precision has been reached by a range of tech¬ 
nical improvements over previous interferometers includ¬ 
ing power stabilization, phase stabilization and increased 
throughput as well as a judicious analysis of detector non- 
linearities, which are the dominant origin of systematic 
error. Furthermore, we have performed the first measure¬ 
ment of third- and fourth-order interference terms, with 
similarly small uncertainties. So far, all our experimental 
results showed no significant higher-order interferences 
and are, therefore, in full accordance with the conven¬ 
tional theory. The dominant sources of imprecision in 
our setup are the uncertainties in determining the de¬ 
tector nonlinearities as well as shot noise in the single 
photon regime. In order to narrow the bound on higher- 
order interference further, highly linear detection systems 
and brighter single photon sources or higher detection ef¬ 
ficiency will be required. A narrower experimental bound 
will aid the development of new theories or constrain free 
parameters of existing ones. In particular, knowledge of 
the various higher-order terms should permit discrimi¬ 
nating between different models for generalized theories, 
such as coefficients in nonlinear extensions of Born’s rule 
[36], the theory of density cubes [37] and quartic quan¬ 
tum theory [38]. For example, the bounds on trans¬ 
late directly to bounds on the magnitude of off-diagonal 
elements in the theory of density cubes [37]. All these 
alternative theories contain quantum theory as a subset 
similarly as quantum theory contains classical theory as 
a subset. The mechanism by which theories exhibiting 
higher-order intereferences reduce to standard quantum 
theory is called hyper-decoherence [38, 39]. This mecha¬ 
nism would be analogous to the process of decoherence, 
which induces the quantum-to-classical transition. Our 
experiment places also a bound on the hyperdecoher¬ 
ence time of the potential extensions of quantum theory 
with second, third and fourth order interference. Such 
post-quantum theories are not only interesting from the 
foundational point of view; they could solve problems 
intractable even on a quantum computer [40]. 
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SUPPLEMENTAL MATERIAL 
SI-ROBUSTNESS OF THE SORKIN EXPERIMENT 

In the following, we will outline the impact of typical error sources occurring in optical multi-path interferometers. 
For the sake of simplicity, we will restrict the analysis to potential second-order interference K 3 arising across three 
paths A, B, C with associated photon transmission rates into the target output mode pA, Pb, Pc and phases 
<j)c- An extension to yet higher orders of interference is straightforwardly possible. There are four types of errors, 
which we are going to consider: 

• Incoherence, e.g. due to phase fluctuation shorter than a single measurement or limited beam overlap on the 
recombining beam splitter 

• Coherent phase fluctuations (phase fluctuations on a time scale longer than each measurement, but shorter than 
a measurement cycle) 

• Input power fluctuations 

• Poissonian photon counting uncertainty (Shot noise). 

It should be reasonable to assume that all sources of errors are independent from one another, such that they can be 
treated separately and their individual influences on the Sorkin experiment can be added. 

We first exclude the counting uncertainty and keep it for later. A general three-path interferometer has re¬ 
spective transmissions Ta,b,c from the input along either of its three paths into the output mode: 

PA = Tapi, Pb = TbP2, Pc = TcPs, 

where the input rates Pi, 2,3 are allowed to be different random variables accounting for power fluctuations between 
the respective single-path measurement settings. 

A two-path measurement of modes A and B takes place at another instance in time and is, thus, subjected to another 
input power P 4 : 


PAB = TaP4 + TbP4 + 2a/T^P4TbP4 cos(^b - <I>a) = Pa{Ta + Tb + 24 /TaTb cos 4)i), (6) 


where the phase difference (j)i = (j)B — <I>A can be treated as a single random variable. Similarly one can express 
the other two-path terms as: pbc = PdiTB + Tc -I- 2\/TbTc cos ^ 2 ) and pac = Pq{Ta + Tq + 2 -y/T^Tc cos with 
independent realizations of phase differences (j)2 = <Pc — 4 >b and <p^ = pA — <Pc (Of course (^1 (j)2 + 4 >s) = 0 must 

hold in a stationary experiment). The two-mode interference term arises as the difference of the two-mode count rate 
and the two corresponding single-path rates: 


Iab = Pab — Pa— Pb 


Pa (Ta +Tb+ 2^/TaTb cos - TaPi - TbP2 


The other two terms are defined accordingly. 


(7) 


The three-path measurement is influenced by yet another input power py and three random realizations of 
the phases 4>a, (Pb, 4>C) which we denote by <(' 4 , 5 , 6 : 


PABC = P7 


Ta+Tb + Tc + 2( \/TaTb cos — (j) 4 ) + a/TbTc cos {pQ — p^) + \/TaTc cos {p 4 — 


( 8 ) 


In this case the three phases influence the measurement at the same time, so their differences can no longer be treated 
as being independent from one another. Of course, the background rate po can also be affected by noise. However, 
as its magnitude is usually much smaller than the other rates, we neglect it in the following analysis. Therefore, the 
second-order interference term es defined in Eq. (4) in the main text is modeled by seven random input powers pi ,,..,7 
and six random phases pi....fi. 

In the experiment we operate at the point of fully constructive interference, that is (pj) = 0 for all phases. In the 
following we will evaluate the impact of the various error sources on the unnormalized term £3 before considering the 
normalization by the sum of the first-order interference terms S 3 . 



Gaussianjfit 



FIG. 6 . Statistical data for a phase fluctuation cr^ = O.Itt. (a) Distribution of the phase differences (/)i, 2 , 3 . (b) Distribution of 
their cosine, (c) Distribution of £ 3 . The red curve is a normal distribution with zero mean value and standard deviation 
for comparison. The ensemble size of the simulation was 10®. 


Phase fluctuations 

We first consider pure phase fluctuations, that is the absence of any changes in input power. Thus pi = pin = const, 
for all 1 = 1,..., 7 and one obtains for the second-order interference term: 

ea = 2pin I \/TaTb [cos (^5 - ^ 4 ) - cos f/fi] -I- \JTbTc [cos (^e - (t>b) - cos (/) 2 ] + V^aTc [cos { 4>4 - - cos 413 ] |. 

(9) 

Whether phase fluctuations influence the experiment in a coherent or an incoherent way will depend on their time 
scale. We will evaluate both cases separately. 


Incoherence 

Incoherence can be modeled by rapid fluctuations of all phases around their mean value {(p) = 0. This means that 
during a measurement the detector integrates over these rapid fluctuations reducing the result from the ideal value 
cosp = 1 to some averaged value {cosp) = X, with 0 < X < 1 measuring the degree of coherence and, thereby, the 
interference visibility. If the incoherence is stationary, it will reduce the interference contrast for all measurement 
settings by the same amount. As evident from Eq. 9, the cosines for each pair of paths enter with opposite signs, 
such that reduced interference contrasts cancel out and no net effect on 63 remains. Thus, the Sorkin experiment is 
immune to such incoherence effects, as has been known for some time [ 1 ]. 


Coherent phase fluctuations 

Slower phase fluctuations preserve the phase within a single measurement, but can alter the phase inbetween mea¬ 
surements of a cycle. This can be modeled by assigning independent random values to the phases (^ 4 , 5,6 adhering to 
a Gaussian distribution with zero mean and a standard deviation of cr,*: (^4 5 e oc An ^2 . The phase differences pi, p2, 
p 3 are differences of two such independent Gaussians. Hence, they are also normally distributed, but with twice the 
variance: 2 3 cx An 20-2 . As in the scenario of incoherence, the phase fluctuations reduce the cosine terms. However, 

’ ’ ' <i> 

not every cosine is reduced in the same way, but by the six independent random phases. One can straightforwardly 
see from Eq. (9) that there is no bias on the unnormalized second-order interference: 

(£3) = 0, 

as all terms are cosines of phase differences with identical spread, so they have identical expectation values which 
cancel each other out. Hence, also coherent phase fluctuations have no systematic influence on our experiment. They 
merely lead to random uncertainties which can be mitigated by averaging over multiple data sets. 

As the cosine is changing its value only in second order in p at the maximum, one can expect a quadratic scaling 
of the uncertainty of €3 with such phase fluctuations. We verify this intuitive expectation by numerical simulations. 
To this end, the random phases pi^...fi are drawn from the Gaussian distributions defined above, with an ensemble 
size chosen large enough to ensure at least two significant digits for all moments in each case. Fig. 6 shows simulated 
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data for a phase uncertainty of = tt/IO in a balanced interferometer {Ta = Tb = Tq = T] p = Tpin as the 
rate transmitted through a single path). As the cosine of a Gaussian variable centered around zero yields a highly 
asymmetric distribution (see subfigure (b)), the resulting distributions for the second-order interference are also 
asymmetric and deviate considerably from a normal distribution. Yet, the expectation value of €3 remains exactly 
zero, as expected from our earlier considerations (c). 

In order to determine, the scaling of the random errors on £3 with the level of phase noise, the procedure is repeated 
for a variety of phase uncertainty levels. The resulting standard deviations are summarized in Table III: For 


50 [tt] 

0.1 

0.03 

0.01 

0.003 

'^^3,ph=,=e b] 

0.69 

0.068 

1 

0 

X 

6.9 X 10"^ 


TABLE III. Scaling of the uncertainty of £3 with phase noise magnitude, 
all cases, there is no bias. The scaling of the random uncertainty can be extracted from the table as being roughly 

C^£3,pha.e oc 7.7pal, 

which is a quadratic scaling, as expected in the vicinity of fully constructive interference. 


Power fluctuations 


Now we consider the other extreme case of perfect phase stability but random fluctuations of the input power. 
Consequently, one can set all cosines in Eqs. ( 6 ) and ( 8 ) to 1 and treat the rates pi ,...,7 as independent Gaussians with 
mean pi^ and standard deviation ap. One obtains for the three-path interference term: 


£3 = Ta {P 7 - P 4 - P 6 + Pi) + Tb {P 7 - P 4 - P 5 + P2) + Tc {P 7 - P 5 - P 6 + Ps) + ■ ■ ■ 


... + 2 


\/TaTb {P 7 — P4) + VTbTc {P 7 — Pb) + \/TaTc {p 7 — P&) 


( 10 ) 


Again, no bias arises as the expectation values cancel each other to zero: 


(£3) = 0. 


Thus, the Sorkin experiment suffers no systematic error from power fluctuations, as long as all shutter settings receive 
the same expected input power. The influence of potential long-term drifts is eliminated by picking the shutter settings 
in random order for each measurement cycle. 

As Eq. (10) is linear in the random powers pi,..., 7 , the second-order interference follows also a Gaussian distribution 
with a standard deviation scaling linearly with Up. The coefficient depends on the specific values of the 

transmittivities. For a balanced interferometer, one obtains 

= 2v^rap = 2v^p^. 

Pin 


Counting Error 

Finally, we take the photon counting noise into account. The number of incident photons in a given time interval 
follows Poissonian statistics, with equal mean photon number and variance. Hence, the count rate measured for 
shutter setting ‘A open’ has a standard deviation of ^/pa, and so on. Again, the second-order interference is unbiased 
by these uncertainties, as the mean values of the individual rates are unaffected. The counting errors for the different 
shutter settings are uncorrelated, so one can add their variances and obtain 

£^€3,count = b/PABC + Pab + Pbc + Pac + Pa + Pb + Pc- 
In case of a balanced interferometer with absolute single path rate p one gets in absence of other error sources 

O'ea.count = 2v^Vp. 

Other sources of noise in the detector can be modeled in the same way and also cause no bias on the result. 
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Combined uncertainty and normalization 


The above discussion shows that none of these error sources has a systematic effect on the unnormalized second-order 
interference. The three sources of random error can be expected to be independent from another, yielding a total 
uncertainty of 

^^3 = + <,powe. + ^"s.count ~ + 11 ((Tp/pi„)") + 2, 

with the numbers in the last step applying to the case of a perfectly balanced interferometer. Evidently, for large 
count rates phase and power fluctuations, which enter in second and first order, respectively, tend to dominate over 
shot noise. Note that the phase noise does not lead to a Gaussian distribution of € 3 , whereas the others do. Hence, 
histogram plots of £3 can reveal which type of error source is dominant. 


It seems more useful to determine higher-order interference relative to first-order interference rather than in 
absolute terms. Therefore, £3 is normalized by the sum of the magnitudes of the first-order terms as defined in Eq.(7): 


_ £3 _ _ 

^3 \Iab\ 


£3 


\Ibc\ + \Iac\ 


( 11 ) 


Clearly, the numerator and the denominator are dependent upon each other. Consequently, fluctuations in any of the 
random variables can lead to an undesired bias on K 3 . 

In order to overcome this problem, we take the averages for £ and 5 separately, as stated in Eq. (5) in the main paper. 
Thereby, we obtain an unbiased estimate for e (which is zero in the conventional theory) and normalize it by the 
averaged first-order interference. The latter is reduced by incoherence and phase fluctuations. The other error sources 
have no systematic influence on 5. 

The relative random uncertainty of the numerator £3 is much larger than the one of the denominator 5^, due to 
(£ 3 ) « 0, but (^ 3 ) > 0. Therefore, the standard deviation of k can be approximated by 


^£3 

(< 53 )' 


( 12 ) 


In case of a perfectly balanced interferometer at fully constructive interference one can evaluate the denominator as 
( 153 ) = 6p, thus: 


i<^p/Pinf) + 2/P- 


S2-ERROR ANALYSIS OF THE DATA 


We conduct a series of M measurements of £ and 6, cycling through all shutter settings. The resulting distribution of 
£ for some exemplary path combinations, normalized by the average of 5, is shown in Fig. 7 as well as in Fig. 3 in the 
main text. As discussed in the previous section, the distribution of the data can provide information on the magnitude 
of random uncertainty as well as on the predominant source of error. In the classical and semi-classical regime the 
central peak of the distributions is slightly more pronounced than expected for an ideal Gaussian distribution. This 
is consistent with the expected shape of a phase-noise induced distribution (cf. Fig. 6(c)). Therefore, the data 
suggests that both, Gaussian power and shot noise as well as non-Gaussian phase noise play a role in these regimes. 
For heralded single photons, on the other hand, the lower count rates render the data shot-noise dominated, causing 
the results to be normally distributed (see bottom row in Fig. 7). 


If the measurement results £^*^ (i = I,...,M) are mutually uncorrelated, one can use the standard error of 
mean to quantify the uncertainty of the average: 

A£=Cre/-\/M, (13) 


with (Te as the measured standard deviation of the data set. The same holds for the uncertainty of n. We verify the 
absence of correlations in the data by calculating their auto-correlation function: 


Aac (fc) = 
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FIG. 7. Histograms for various orders of interference and measurement regimes. Top row: Data obtained in the classical regime 
with M = 5618 measnrements for paths ABCD (a) and ABCDE (b). The corresponding three-path interference data is 
shown as Fig. 3 in the main paper. Middle row: Data obtained in the unheralded single-photon regime (M = 1912) for paths 
ABC (c) ABCD (d) and ABCDE (e). Bottom row: Data obtained in the heralded single-photon regime (M = 1912) for 
paths ABC (f) ABCD (g) and ABCDE (h). The blue lines are Gaussian hts of the expected distributions of the data and 
the black bars indicate the mean values with their widths as the standard error of mean. 


An exemplary auto-correlation for £3 measured with the laser source on paths A,B,C is shown in Fig. 8(a). 
Evidently, no significant correlation among the data points can be detected, not even for subsequent ones. The 
same observation is made for all other measurement settings and light sources. Therefore, it is fully justified to 
use the standard error of mean, as defined in Eq. (13), as a measure for the uncertainty of the higher-order interference. 

One has to be more careful, however, when it comes to averaging across the various path combinations in 
the multi-path interferometer. For example, phase fluctuations in a single path will influence the interference 
among several path combinations simultaneously. Therefore, cross-correlations between the measured data for 
different path combinations {i = ABC ,..., CDE for three paths) can be expected. Indeed, one finds non-vanishing 
cross-correlations {Rcc)ij for £3 and £ 4 , as shown in Fig. 8(b) and (c), respectively. With this at hand, one can 
obtain the standard error of mean of the final result for K 3 via error propagation: 


Ak3 


10 ^ 


CDE 

E AK3yA«;3j , 

J^ABC 
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FIG. 8 . Auto- and crosscorrelation of the experimental data obtained with the laser source, (a) Autocorrelation of ea for path 
combination ABC. Only the first quarter of the correlation function is shown. For nncorrelated normally distributed data 
95% of the points should lie inbetween the dashed lines, (b) Crosscorrelation matrix between the second-order interference ea 
measured for the various path triplets, (c) Same for the third-order interference £ 4 . 


with denoting the standard error of mean of K 3 on path combination i. Analogously, one gets for K 4 : 


Ak 4 = 




BCDE 

E Ak 4 iAK4 j. 

i,j=ABCD 


These are the experimental uncertainties presented in Table I in the main text. 


S3-TOMOGRAPHY DATA AND PREDICTED COUNT RATES 


For completeness we show additional data/figures resulting from our measurement: 

The measured count rates of the different path combinations for the measurement with the laser can be seen in 
Fig. 10. After filtering the 5618 measurement sets recorded in a measurement time of 68 h the standard deviation 
was measured to be 0.3%. A bar chart with the average powers over the whole measurement duration can be seen in 
Fig. 2 of the main text. 

The measured count rates of the different path combinations for the measurement with the single photon source 
for the semi-classical measurement can be seen in Fig. 10. After hltering the 1912 measurement sets recorded in a 
measurement time of 88 h the standard deviation was measured to be 3.6%. This higher standard deviation compared 
to the classical light source results mainly from shot noise and due to the fact that the power of the blue pump laser 
was not stabilized. 

Fig. 11 shows the comparison between the measured values (shown in Fig. 2 in the main text for the laser) and the 
expected values arising from the reconstructed density matrix for all different path combinations. 


S4-NONLINEARITY MEASUREMENT AND ASSOCIATED UNCERTAINTY 


Detector nonlinearities are the predominant cause of systematic deviations in our implementation of the Sorkin 
experiment. We deduce these nonlinearities from independent beam combination experiments [2, 3]. Their expected 
impact on the higher-order interference term nth is then calculated by applying the nonlinear detector response function 
to the various multi-path probabilities calculated from the density matrix of the interferometer. In the following, it 
will be discussed how and with which uncertainty the nonlinearities are reconstructed and how this translates into 
the theoretical prediction Kth and the error bars of this prediction. 
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FIG. 9. Measured intensities for the different path-combinations in the classical regime. Each combination is indicated by a 
separate color within each panel. 
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FIG. 10. Measured intensities for the different path-combinations in the semi classical regime. Each combination is indicated 
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FIG. 11. Comparison between measured data and data simulated from tomography for (a) the classical and (b) semi-classical 
measurement. 


Nonlinearity of the photodiode and the voltmeter 

Polynomial expansion of the transfer function 

In a measurement instrument for classical light the output voltages V are related to the impinging optical powers 
p by some transfer function p = f{V). A variety of effects can contribute to nonlinearities in a photodiode and the 
attached voltmeter, such that it is difficult to develop a specific model. However, as the nonlinearities are quite weak 
compared to the linear part of the transfer function (on the 10 ppm-level), one can reasonably expand f{V) into a 
polynomial of degree n in the region of interest [3]: 


i=o 

with coefficients a = (oq, ..., a„)^ and the residual r„(I/, a) quantifying the fidelity of this approximation for a given 
measured voltage V and polynomial degree n. 

The goal of the beam-combination method is to find a suitable polynomial for f{V). One overlaps two inco¬ 
herent beams, which can be shuttered and power-controlled individually, on a detector (see, e.g.. Fig. 1 in [3] for an 
illustration). For a series of k = 1,... ,M settings of optical powers, the voltages Vi^k, 1^2,fe for either of the paths 
being open, Vs^k for both path being open and Vb.fc for both paths being closed are recorded. The corresponding 
unknown optical powers are Po.fcj ■ ■ • iPa.fc- For incoherent beams, the optical powers should be additive, i.e., 

yk : p3^k + Po,k = Pi,k + P2,k- (14) 

Condition (14) leads to the following constraints on the coefficients: 

n 

Vfc : ^ a, (Vl, + VI, - VI, - = r„(Vk, a), (15) 

with r„(Vfe,a) = r„(Vi^fc,a) +'r„(V 2 ,fc, a) — r„(V 3 fe,a) — r„(Vo,fej3-) measuring the combined residual for data set k. 
Note that the coefficient oq drops out of the sum due to Vq, = ... = V^, = 1. This is to be expected, as constant 
off-sets should have no influence in a balanced zero-sum experiment as given by Eq. (14). 

The aim is to determine the nonlinearity of the detector. The exact value of the linear slope oi is of no 
interest in this respect. Thus, it can be set to Oi = 1 for simplicity. It is convenient to introduce the quantity 

S'j-fe = vl, + vl, - vi, - vl,. 











15 


With this at hand, one can reformulate the constraints (15) into the following least-square optimization problem: 

Rn{a) = ^ J -2-; Rnia) Min, (16) 

k=i 

which is solved by finding a configuration of nonlinearity parameters a = ( 02 ,. •. ,a„)^ that minimize i?„(a). Here, 
the contribution of each dataset k is inversely weighted with its variance These variances are a result of random 
experimental uncertainties of the measured voltages and may vary across the measuring range of the instrument. 
Taking them into account bears the advantage that more certain parts of the data contribute more to the optimization 
than uncertain ones. 

If these data uncertainties are known, one can straightforwardly solve the least-square problem [4]: 

a = CA^b, 

with the covariance matrix 

C=(A^A)~', (17) 

the weighted matrix elements = Sj^k/<^k (for k = 1,.. . ,M and j = 2,... ,n) and vector components b^. = 

—Si^klo'k- Then, the optimal transfer function reads 


fopt{V) = V + Y,ajV^, (18) 

i=2 

where the irrelevant off-set term ag has been set to zero for convenience. 


Uncertainty of the nonlinearity 


In the following, it will be explained how the variances cr^ can be extracted from the data and how they influence the 
uncertainty of the optimization. The experimental data of [3] is used for this procedure and the subsequent prediction 
of Kth- 

The optical power in the beam combination experiment is increased linearly between subsequent measurement points 
(see Fig. 12(a)). The uncertainty of each voltage Idn.fc is estimated by calculating a floating standard deviation cry^ ^ 
of the surrounding data points (r = 101 points were taken in each case) with respect to a perfectly linear slope 
(Fig. 12(b)). The uncertainty ak can then be taken as the standard deviation of each summand in Eq. (16), i.e., 
as the standard deviation of the function y(Vk,a) = Si^k + ^ 2^=2 4 voltages in each dataset Vk are 

independent from one another and if the uncertainties of the voltages are small compared to their mean values, one 
can express as [5]: 



ay(Vk,a) 



i=2 



In general, this depends on the measured voltages as well as on the nonlinearity coefficients. However, for a weak 
nonlinearity ajV^ <C V holds for all j > 1. Hence, ajV^~^ <C 1 and ak is dominated by uncertainties in the linear 
part of the optimization problem: 


m—0 

The resulting uncertainties are shown in Fig. 12(c) and are used as weights in the least-square optimization (16). 

The variances and covariances of the resulting coefficients a can then be readily obtained from the covariance 
matrix (17) [4]: 

^aj ~ — — Covi^dj J = Cj-i^k—l- 
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FIG. 12. (a) Measured voltages in the beam combination experiment, (b) Uncertainties of these voltages, obtained from a 
floating standard deviation, (c) Total uncertainty of the dataset, which is used as weight in the least-square optimization. 
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FIG. 13. (a) Normalized residual X{n) of the least square optimization in dependence of the chosen order of the polynomial 
n. (b) Resulting nonlinearity and associated uncertainty for n = 3. 


The uncertainty of the reconstructed transfer function in Eq. (18) at a given voltage V can be calculated from the 
covariance matrix via error propagation 




(r)= E 

p,q=2 


dfopt dfopt 
dur) dua 




p,q=2 


(19) 


Note that this purely data-based approach makes no assumption on the precision of the instruments or the magnitude 
of random fluctuations in the experimental conditions. 


Order of the polynomial and resulting transfer function 

The optimization procedure Eq. (16) is solved for various degrees n of the polynomial. We select the lowest possible n 
beyond which the quality of the optimization does not significantly increase any further. To be specific, we calculate 
a test quantity X{n), defined as the residual i?„(a), normalized by the remaining number of degrees of freedom in 
the system X{n) = Rn{a)/{M — (r — 1) — n -b 1), with M — {r — 1) being the number of data points remaining after 
calculation of the local standard deviation. 

This quantity has been calculated for increasing n and the result is plotted in Fig. 13(a). As one would expect, a 
linear function n = 1 fares worse than nonlinear ones, which suggests that indeed significant nonlinearity is present in 
the system. There is some improvement in the quality of the ht from n = 2 to n = 3 {X{i)/X{2) « 0.98). However, 
beyond n = 3, the optimization improves barely any further: A'(4)/A(3) « 0.992 and A(5)/A(4) « 0.998. This 
suggests that n = 3 is a reasonable choice for the order of the polynomial without overinterpreting the data. This 
value is chosen for the estimation of the nonlinearity in this work. 

The resulting parameters, which optimize (16) for the beam combination measurements of Ref. [3] are: 


a = (-59.8,3.1)^ X 10"®, 


(20) 




















with the covariance matrix 
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f 496 -45 
-45 4.2 


X 10"^^. 


( 21 ) 


Implications on the Sorkin experiment 

Now it will be discussed how the systematic influence of nonlinearity-affected detection of the interferometer’s output 
powers on higher-order interferences can be quantified. To this end, the state of the interferometer, i.e., all amplitudes 
and phases of the 5 paths, is reconstructed via state tomography. From this the expected optical powers of the Sorkin 
terms {pa,Pab, • • ■) are calculated (see Fig. 11). Due to the nonlinear dependence between these powers and their 
measured voltages (Va = foptiPA), Vab = foptiPAB), • ■ •) an apparent non-zero higher order interference arises, 
even if no physical higher-order interference is present. Its value will be predicted by calculating the nonlinearity 
affected voltages via the inverse transfer function. 

Clearly, any uncertainty in fojtip) translate directly into an error bar of the predicted Kth- The inverse function 
can be analytically calculated from fopt(V) for low-order polynomials. However, the measurement device is mostly 
linear, such that p = f{V) = V + /nl(C), with |/nl(C)| <C V. With this one can approximate 

/■^(p) /nl(p) = 2p-/(p). (22) 

Then, clearly, the uncertainty of the inverse function is approximately equal to the one of the forward transfer function. 

af~i{p) « CT/(p), 

which can be calculated via (19). The inverse transfer function and its uncertainty are shown in Fig. 13(b). 

Note that the tomography data itself will also be affected by the nonlinearity, so one could apply the func¬ 
tion fopt{V) to reconstruct the ‘true ’optical powers before the density matrix is reconstructed. However, the 
higher-order interference must always be zero for any physical density matrix, regardless of its precise structure (if 
Born’s rule holds). Therefore, a forward nonlinearity correction of the tomography data, leading only to minuscule 
changes in the density matrix, has no significant impact on Kth and can, hence, be omitted. 

It is of course difficult to tell in which way the realizations of f~^{p) for the various values of p are corre¬ 
lated with each other. One may consider two extreme cases: The case of ‘maximum correlation’, where all values of 
foptip) are shifted in the same direction by one standard deviation and the completely uncorrelated case, where each 
value fojtip) is replaced by an independent Gaussian random variable with mean fop!tip) and standard deviation 
(T/(p). In the maximally correlated scenario, one finds generally much smaller error intervals than in the uncorrelated 
case. Therefore, we resort here to the latter scenario to provide the more conservative error estimation. 

For such a fully independent random variation of the nonlinearity for each power term, one can calculate the uncer¬ 
tainty of the prediction of Kth by propagation of Gaussian errors. With this procedure one obtains for the dataset 
of the main paper and the corresponding tomography data, the apparent higher order interferences shown as the red 
symbols in Fig. 14. Similar values and uncertainties are expected for all subsets of the five interferometer-paths. In 
case of second- and third-order interference, one can average over these subsets to obtain, along with the prediction 
for the fourth-order interference on the full set of paths, the final estimates used in the main text of the paper: 

= (9.7 ±3.1)-10-5 
4th = (-1-6 ±4.1)-10-5 
4 th = (-3-9±5.1)-10-^ 

The uncertainties of the measured value {k) and the predicted value Kth arise from independent experiments. There¬ 
fore, it is natural to assume that there is no correlation between them, which implies that the uncertainty of the 
corrected higher-order interference k = (k) — Kth is obtained by adding the individual variances: 

Ak = Ak^ + 

These are the error bars shown in Table H and Fig. 5 in the main text. 
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FIG. 14. Experimentally measured higher-order interferences (blue data points and error bars) in comparison to the predicted 
values Kth arising from detector nonlinearities (red). All possible path combinations are shown. Errors bars indicate one 
standard deviation, (a) Second-order interference K 3 . (b) Third-order interference K 4 . 


Nonlinearity of a single photon detector 


The reverse-biased avalanche diodes used as single photon detectors have an intrinsic deadtime t. The resulting 
saturation of the detector for increasing count rates is the dominant nonlinear mechanism. The transfer function 
p = f(V) (p and V now being count rates) and its inverse can be expressed by a simple saturation model [3]: 


f(V) 


V 


1 -tV 


; r\p) 


1 

1 + rp 


(23) 


Just as for classical light detectors, the unknown parameter r can be determined via beam combination experiments. 
The additivity of the optical powers (14) translates to the following condition: 

Vfc : F{t, Vk) =/-1 [/ {Vo,k) + f [/ + f{V2,k)] 

_Vo,k + Vs^k — 2 TVo,fcV 3 _fc bi.fe + b2.fe — 2TVi_fcV2,fc _ 

1 — 1 — T'^Vi^kV2,k 

The optimal solution for the deadtime is found via the nonlinear optimization problem: 


R{t) = Y, ^ Min. (24) 

fc=i 

As before, is the variance of each function F(r, Vk) and is used to weight the optimization. If one is sufficiently 
far away from the saturation point, the nonlinearity is weak and the variance is dominated by the uncertainty in 
the linear part of F(t, Vk). Then, cr^ « classical detector. As before, the individual 

uncertainties cry^ ^ can be calculated from the data via a floating standard deviation. 


Using this technique, we obtain a deadtime r = 33.9 ns with an uncertainty ar = 0.3 ns for the detectors 
used in our experiment. One can then use the saturation model (23) to calculate a prediction of higher-order 
interference from the tomography data, just as for the classical light measurements. In this case, however, the 
uncertainty in r must influence all count rates in the same direction. Therefore, the uncertainty in the prediction can 
be estimated by performing it also for deadtimes tFut- These predictions are shown together with the experimental 
data in Fig. 15 for the individual path combinations. Compared to the classical case, they have negligible error bars, 
as a very specific nonlinearity model with a single, precisely determined parameter is used. The averages over all 
path combinations yield: 


= (-11.18 ±0.10) X 10"^ 
= (-3.48 ± 0.03) X lO"'^ 
4% = (2.85 ±0.05) X 10"®. 
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FIG. 15. Same as Fig. 14, but for the single photon data. Error bars of the prediction are invisible on this scale. 


Nonlinearity for the heralded single photons 


In case of heralded single photon detection one has to take into account that both, the heralding detector (h) as well 
as the detector fed by the interferometer (i) are inactive for the dead time r whenever they detect a photon. Some 
of these cases are coincidences (c), which blind both detectors. One can derive a saturation model, with the ‘true’ 
coincidence rate Pc depending on the measured rates Id.h.c- 


Pc = f{Vi,Vi„V,) 


K 

i-T{Vi+Vi:-v,)' 


This model is applied to correct the coincidence data obtained in the heralded single photon regime, which is then 
used to directly calculate the corrected higher order interference presented in Eq. (11) in the main text. As 

for the unheralded single photons, the error interval of the correction is estimated by performing it with deadtimes 
T ± (Tt-. Again, this correction uncertainty is assumed to be independent from the experimental errors in the Sorkin 
measurement, such that their variances can be added. 
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